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AN ANALYTICAL METHOD FOR DETERMINING 


THE MELT-LAYER THICKNESS AND REGRESSION RATES 
FOR A HYBRID OR TRIBRID ROCKET MOTOR* 

By Walter Brooks Drew 
Langley Research Center 

SUMMARY 

An analytical investigation was performed that produced a method whereby the 
regression rates and melt -layer thickness of a tribrid solid fuel could be determined. 
Slab-burner tests of the fuels in question substantiate the analytical method. Regression- 
rate studies of hybrid fuels in the past have not been concerned with the heat transferred 
into the grain because previous solid fuels have had low thermal conductivity and did not 
allow appreciable heat penetration into the grain. With the appearance of hybrid and 
tribrid fuels that contain a large amount of high conductivity metal, a need arose for a 
method of predicting regression rates that took into account a new problem that pre- 
sented itself, melt-layer buildup. In order to predict melt-layer buildup, the true nature 
of heat transfer, liberation, and absorption in the fuel grain must be considered. The 
method presented has accomplished this objective. 

The investigation was originally undertaken to demonstrate the effect that an obvi- 
ously important parameter, fuel thermal conductivity, has on melt-layer buildup and 
regression rate. Also, an interesting result of this analysis was the delineation of the 
effect that heat input to the fuel surface has on melt -layer buildup. 

Analytical results based on laboratory-size motors indicate that the slab-burner 
and cylindrical -grain configurations will produce melt-layer thicknesses greater than the 
experimentally determined critical value for stripping of 0.125 inch (0.32 cm) and will 
permit the liquid fuel to be stripped away from the fuel surface by the viscous forces in 
the gaseous boundary layer. This condition will produce low combustion efficiencies in 
the tribrid rocket system. 

Calculations also show that the thickness of the melt layer is a strong function of 
the convective and radiative heat input to the fuel surface and can be minimized by manip- 
lilating surface heat input through flow rates and grain geometry. 

*The material presented herein was included in a thesis submitted in partial fulfill- 
ment of the requirements for the degree of Master of Mechanical Engineering, University 
of South Carolina, Columbia, South Carolina, 1968. 



INTRODUCTION 


The hybrid rocket propulsion system utilizes a solid-fuel grain with a liquid oxi- 
dizer sprayed across the grain surface. Research has developed the technology neces- 
sary for the production of a flight hybrid system, and the hybrid concept has been satis- 
factorily demonstrated. 

The "tribrid" rocket propulsion system is a natural extension of the hybrid tech- 
nology. The difference between the two concepts lies in mixing a third nonworking gas 
component into the hot hybrid combustion products. The best fluid to be added in this 
manner is hydrogen, because of its low molecular weight. The consequences of adding 
hydrogen to the hybrid system are advantageous in two ways. First, the temperature 
of the gases entering the rocket nozzle is reduced from above 8000° F (4426° C) to 
below 3000° F (1649° C) for this study. The advantages of lower nozzle temperature 
appear in simpler fabrication techniques, lower weight, and better reliability. Second, 
the specific impulse is increased since it is proportional to the square root of the stag- 
nation temperature divided by the molecular weight of the exhaust gases, and at some 
percentage H2, these parameters will produce a maximum specific impulse. Figure 1 
shows a schematic of a tribrid rocket motor. 

The particular tribrid system under study in this investigation is composed of 
liquid fluorine as the oxidizer, lithium as the fuel, and hydrogen as the nonworking com- 
ponent. The theoretical specific impulse Is for this system is 511 seconds, as com- 
pared with 478 seconds for the liquid hydrogen, fluorine rocket which represents the 
highest yield liquid rocket motor, and 440 seconds for high-performance hybrid motors. 
This increased specific impulse is expected to more than offset the additional hydrogen 
tankage requirements. 

Preliminary slab burner tests of the tribrid system performed by R. J. Muzzy of 
United Technology Center indicate that a major problem exists when employing a lithium 
fuel as the solid grain. This problem occurs when the lithium fuel melts faster than it 
can be vaporized and burned. After a 0.33-cm-thick layer of melted lithium forms on the 
surface, it will be stripped away by the viscous forces of the boundary -layer gases. When 
this condition occurs, the unburned lithium is exhausted and high inefficiencies occur in 
the system. This liquid layer is produced on the lithium fuel because the fuel has a high 
thermal conductivity at a high temperature (-2400° F (1315° C)). Also, lithium has a 
large heat of vaporization as compared with the heat of fusion (186 Btu/lb (103.3 cal/g)) 
at a temperature of 356° F (166° C). These physical characteristics allow energy to be 
transferred into the grain at a rate larger than the rate at which the surface is burned 
away. 
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A proposed solution to this problem is to decrease the thermal conductivity of the 
fuel grain by mixing the lithium with a low thermal conductivity hydrocarbon solid fuel 
and adding an appropriate amount of oxygen to the fluorine in order to burn this fuel. 

This approach to the problem appears to be feasible, but it has an inherent drawback 
because of the fact that the more hydrocarbon binder that is added, the lower the specific 
impulse will be. Also, after a certain amount of binder has been added, self-sustaining 
solid reactions will occur. Anal 3 d;ical studies have shown that a 90 -percent lithium, 
10-percent binder is theoretically capable of 502 seconds; an 80-percent lithium, 
20-percent binder will produce 491 seconds specific impulse; and if the lithium content 
is reduced to 70 percent, the specific impulse will drop to 480 seconds.* From these 
numbers the importance of determining the minimum amount of binder required to pre- 
vent melt-layer buildup is appreciated. 

Herein lies an important objective of this study — to determine analytically the 
minimum amovmt of binder necessary in order to prevent the formation of a melt layer 
on the tribrid fuel grain surface that could be stripped away by the viscous shear forces 
of the reacting boundary layer. 

It is obvious that the melt -layer thickness is directly related to the difference 
between surface and melt-layer regression rates. Regression-rate studies of hybrid 
fuels in the past have not been properly concerned with the internal heat transferred into 
the grain because previous hybrid fuels have had low thermal conductivity and did not 
allow an appreciable amount of Internal heat transfer by conduction. In order to predict 
melt -layer buildup on a lithium fuel, all aspects of heat transfer in the fuel grain must 
be considered. 

An analytical solution to this problem, which is dependent upon experimental values 
for surface temperature, is presented. The solution was accomplished by using estab- 
lished methods to determine the heat transferred from the turbulent reacting boundary 
layer by convection and radiation to the fuel surface. The analysis considers the heat of 
vaporization of the fuel at the surface, the conduction through a melt layer to the solid- 
melt interface (the heat of fusion being taken into accoimt at this point) and finally, con- 
duction into the solid grain. 

This investigation was originally imdertaken to demonstrate the effect that an obvi- 
ously important parameter, fuel thermal conductivity, has on melt-layer buildup and 
regression rates. The results of the analysis also made clear the influence of the heat 
input to the surface from boundary-layer convection and radiation. Details of the method 
are presented herein. 


*This information was obtained informally from impublished data supplied by the 
Propulsion Branch of the Langley Research Center. 
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SYMBOLS 


A area, in^ (cm^) 


Ai = At k (--^ + dimensionless 

" \ Ar2 2ri Arj’ 

B mass transfer number, dimensionless, ^ ^ 

Ub % 

Bi constant, 1 + dimensionless 

Ar^a 

B’ = B X dimensionless 


C constant 


Cf skin-friction coefficient, dimensionless 


D 

known constants of simultaneous equations 


H 

sensible enthalpy, Btu/lb (cal/mole) 


AH 

sensible enthalpy difference between flame and wall. 

Btu/lb (cal/mole) 

K 

mass fraction of nonvaporizing matter in^/ sec 

PCp 

(cm2/sec) 

2 At K 
" irk 

, sec-°F-in2/Btu (sec-°C-cm2/cal) 


Ko = 1 + ^ dimensionless 

Ar^ 


2 Atk 
Ks = 9 , 

Ar^ 

dimensionless 


Ko,e 

oxidizer concentration in free stream 


Ly 

heat of vaporization of vaporizing component, Btu/lb 

(cal/mole) 

N 

radiation parameter, dimensionless; also as integer 


Npr 

Prandtl number, dimensionless, 

k 
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NRe 

Nst 

O/F 

Q 

R 

^oj^bs 

T 

Z 

Cp 

dm 

h 

hv 

k 

m 

r,e,z 

r 

t 

u 

x,y 

a 


Rejniolds n\unber, dimensionless, 

Stanton number, dimensionless, ^ ^ „ 
oxidizer to fuel ratio 

heat transfer, Btu/in^-sec (cal/cm^-sec) 
radiation, Btu/in2-sec (cal/cm^-sec) 
initial and back side radii, respectively 
temperature, (°c) 
radiative path length, in. (cm) 
specific heat, Btu/lb-°F (cal/deg mole) 
melt-layer thickness, in. (cm) 

convective heat -transfer coefficient, Btu/in^-^F-sec (cal/cm^-^C-sec) 

heat of gasification of vaporizing component, Btu/lb (cal/mole) 

coefficient of thermal conductivity, Btu/in-sec-°F (cal/cm-°C-sec) 

mass-transfer rate, lb/sec-in2 (g/sec-cm^) 

cylindrical coordinates 

regression rate, in./sec (cm/sec) 

time, sec 

velocity, in./sec (cm/sec) (in x direction) 
rectangular coordinates, in. (cha) 
empirical radiation constant, dimensionless 
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emissivity, Btu/in^-sec (cal/cm2-sec) 

(I viscosity, Ib-sec/in^ (g-sec/cm^) 

p density, Ib/in^ (g/cm^) 

Pjj density of nonvaporizing component, Ib/in^ (g/cm3) 

a Stefan -Boltzmann constant, 0.1713 x 10"® Btu/ft2-hr-°R^ (cal/cm2-sec-° C4) 

T skin friction, lb/in2 ^g/cm2) 

Subscripts: 

b flame zone 

bs back surface 

c convection 

cs conducted from siirface 

e edge of boundary layer 

f fuel 

in radiation and convection to the surface 

I liquid 

m melt 

o initial value 

r radiation 

Ar space increment 
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s 


solid 


V vaporization 

w wall 

X x-direction 

t time 

At time increment 


ANALYSIS 

The basic theory involved in this study can be divided into two parts: First, the 
heat convected and radiated to the fuel surface must be established by considering a 
turbulent reacting gaseous boxmdary layer above the solid fuel. Second, the utilization 
of this energy in the fuel grain is treated in detail, generally through the use of numerical 
methods. Means for determining the physical parameters associated with the gaseous 
boundary layer and solid fuel are given. A computer program was produced that will 
make the analysis more readily usable. 

Method for Calculating Heat Input to Fuel Surface 

Marxman and Wooldridge have shown transition between laminar and turbulent flow 
occurring at = 10^ (refs. 1 and 2) for the typical hybrid application. This transi- 

tion produces turbulent flow under typical conditions. They assumed steady-state condi- 
tions where all the heat transferred to the surface is balanced by the mass -flow rate of 
fuel coming from the surface multiplied by the total heat of gasification of the fuel; that 
is. 


Qw — rnf^v ~~ 


( 1 ) 


where 

nif 

hv 

Pf 


mass rate of fuel leaving surface 
an effective heat of gasification 
density of fuel 

regression rate of fuel surface. 
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Although hy is composed of the heat of vaporization and the heat required to 
bring the internal energy of the grain up to the vaporization temperature, the true nature 
of heat conduction into the grain has been avoided. (See refs. 1 and 2.) In past studies, 
this assumption has been justified by the fact that solid-fuel rocket grains are thermal 
insulators. 

The heat into the wall is divided into two parts: that due to convection and that due 
to radiation. Convective heat transfer in a turbulent reacting boundary layer is discussed 
first. In order to clarify the terminology initially, let us consider figure 2. 

Convective heat transferred into the grain surface can be expressed as 


^CW Cp\9yyl 


W 


( 2 ) 


In equation (2), it appears that the dimension y has been normalized. At this point the 
Stanton number Ngt is defined as 


Not = — — 
Vb“b 


and equation (2) becomes 


^cw - NstPb% 


(3) 


where Ngt is based on a AH that is equal to the sensible enthalpy difference between 
the flame and the wall ^AH = (cpT^j^ - (cpT)^, Btu/lb, where Cp is the specific heat of 

the gases at the point in question^. 

Combining equations (1) and (3) yields an expression for the regression rate 




(4) 


Equation (4) can be made usable by expressing Ngt in terms of the skin-friction coeffi- 
cient Cf since experimental data and empirical expressions are available for Cj. To 
do so, the Reynolds analogy between heat transfer and the local skin friction is used. 

For a nonunity Prandtl number, this relationship takes the modified form 


Qw 

AH 


= — Npr 
Ub 


-0.67 


( 5 ) 
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From the definition of Cj, 


Cf = 


2t, 


w 


PpU. 


e“e 


Substituting equation (5) into equation (6) yields 


^PeUe2 
2 

and equation (7) into equation (3) yields 


. Cf PgUe^ -0.67 * „ 
Qw=Y-sg-Npr ah 


■KT « „ A TT ^f ^e^e XT -0.67 * „ 
NstPb% AH = Y" ' ub ~ ^Pr AH 


Now an expression for Ngj in terms of Cf/2 can be formed, 

2 


XT ^f Pe^e XT 

Nst = -o o Npr 

^ PbV 


-0.67 


( 6 ) 


(7) 


( 8 ) 


(9) 


Substituting this expression for Ngj- into the regression rate (eq. (4)) 3delds 


r 


^f Pe^e^ 
^ Pb%^ 


Npr 


-0.67 Pb ^ 
Pf b hy 


( 10 ) 


At this point it is convenient to define a mass transfer munber B which has been modi- 
fied to account for the effects of a reacting boimdary layer (ref. 3). 


B - ah 
Ub hy 

and the regression rate equation (eq. (10)) becomes 


PeUe Cf -0.67 
Pf 2 


B 


( 11 ) 


Because the definition of B takes into account physical properties of both the 
solid and gas phase, it is regarded as a thermochemical constant that characterizes the 
particular propellant. And further, B' is defined as B x Npr" * and this param- 
eter is the similarity parameter of a boundary layer with mass injection. (See ref. 2.) 
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Because the Prandtl number of a gas is near unity for a turbulent boundary layer, B 
serves a dual role as a thermochemical and boundary-layer flow parameter. In order 
to put equation (11) into a useful form, the friction coefficient Cf must be determined 
for a tiirbulent boimdary layer with mass addition and combustion. 

The turbulent skin-friction coefficient for a Newtonian fluid is given by Schlichting 
(ref. 4) as 

^ = 0.03NRe_x‘“'^ (12) 


Marxman and Wooldridge (ref. 2) have shown that the effects of mass injection upon the 
friction coefficient can be expressed by 


Cf,o L 


loge(l + B-) 


nO.8 


B’ 


1 + 


13B« ^ 
10 


11 


(1 + BO (1 + |- 


i\2 


(13) 


and a much simpler formula can be used in the range 5 < B’ < 100 which tjq)ifies hybrid 
operation 


Cf 



1.2B’ 


-0.77 


(14) 


Combustion affects the skin-friction coefficient by creating variable properties in 
the boundary layer. Investigators (refs. 5, 6, 7, and 8) have shown that this effect can be 
accounted for by correcting pg by a density ratio of the form (p/Pe)^’^ where p/Pg 
is a suitably defined reference density. For most hybrids this ratio can be taken as 
unity. (See ref. 8.) 

By substituting equations (12) and (14) into equation (11), the final expression for 
the regression rate is obtained: 

i 0.036/' p\®*®„ „,0.23 ^t -0.2 

r = NRg^x (15 

Two basic assumptions have been made in deriving equation (15): 

(1) The applicability of the Reynolds analogy 

(2) That a boundary-layer treatment is valid even for the relatively high rates of 
surface mass flux encoimtered in hybrid operation. 
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Assumption (1) was justified by Wooldridge and Muzzy in reference 9. These 
investigators showed that the turbulent Prandtl number, Schmidt number, and therefore, 
the Lewis number were unity. This statement was verified both anal 3 dically and 
experimentally. 

As far as assumption (2) is concerned, velocity profile measurements by Mickley 
and Davis (ref. 10) and Wooldridge and Muzzy (ref. 9) show excellent agreement with the 
semiempirical, explicit -form velocity-profile equations up to values of B' of 100. All 
practical hybrid or tribrid rocket motors fall within this range. 

At the present time there are many controversial aspects concerning the theoreti- 
cal model and the physical reality, but it is the opinion of the author that the presented 
approach is the best available. The radiative heat transfer to the grain smrface has been 
described by Marxman (refs. 2 and 8) as 


Qr = ae^T/(l - 


( 16 ) 


where 

o Stefan-Boltzmann constant 

€Tff emissivity of wall 

Tj, effective radiation temperature 

a empirical radiation coefficient 

N radiation parameter 

Z radiative path length 


The form of equation (16) has had recent modifications; however, the presented form is 
considered to be adequate for this particular application. 

The energy reaching the fuel surface cannot be expressed by a simple addition of 
the heat due to convection and radiation. A complication arises because as the surface 
absorbs the radiative heat, this heat increases the regression rate and therefore the 
blocking effects tend to reduce the convective heating. The following formula has been 
shown (ref. 8) to express this relationship 


Pf 




(17) 
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where is the convective heat transfer in the absence of radiation. It might be 

noted that for small values of trade-off is nearly exact and the regression 

rate may be calculated with little error by using alone. 

Many hybrid and tribrid rocket systems contain a high percentage of metal loading 
in the grain. This metal can create significant radiative heat transfer if its combustion 
products are solid, but if the combustion products are gaseous, little radiation will be 
exhibited; therefore, the effects of metal particle loading have been investigated for both 
radiative and nonradiative systems. 

Reference 3 indicates that mass addition of particles in the flow has little effect 
on the shear stress distribution or velocity profile shape inside the boundary layer. This 
statement means that the addition of particles from the surface will not affect the relation- 
ships between Cf and Cf^Q. The blocking effects which reduce Cf should be depend- 
ent only on the gas phase inside the boimdary layer as long as the solid particles occupy a 
negligible volume compared with that of the gas. This condition would have to be satis- 
fied in any real situation. Because of this requirement, the regression-rate equation 
should be written in terms of the component in the grain that produces these boundary- 
layer gases. This phenomenon can be taken into account by defining an effective heat of 
gasification as (ref. 3) 


PvK,eff = Pv^^r,b + ^m (18) 

where 

Py density of vaporizing component 

Cjjj specific heat of nonvaporizing component 

AT temperature difference between surface and ambient deep in grain 

K mass fraction of nonvaporizing matter 

hvjb gasification of vaporizing component 

Equation (18) illustrates the fact that one component is heated to the surface temperature 
and vaporized, whereas the other component is heated only to the surface temperature. 

Generation of a Regression-Rate Temperature-Profile Matrix 

The preceding section gave a method whereby the heat input to the grain surface 
could be calculated. Now the problem is to represent mathematically the distribution of 
this heat energy inside the fuel grain. The physical situation is as follows: 
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(1) Fuel is vaporized at the surface by absorbing that amount of heat of 
vaporization. 

(2) Energy is transferred across the melt layer and absorbed by increasing the 
internal energy of the fuel in the melt layer. 

(3) At the melt-solid interface the heat of fusion of the fuel is subtracted from the 
total heat input from the surface. 

(4) Energy is transferred by conduction from the melt-solid interface into the 
solid grain and thus the internal energy of the solid fuel is increased. 

Carslaw and Jaeger (ref. 11) have stated that apart from the few exact solutions 
(concerning phase change) all problems have to be attacked by numerical methods, and 
that numerical techniques have the advantage that the variation in the thermal properties 
with temperature can be taken into accoimt. A numerical solution was utilized for this 
analysis and it proceeds as follows: 

The partial differential eqioation that describes the heat flux through a cylinder with 
d^T/de^ =0 is 

0T _ 1 3T 

at - r 8r q^2J 

Experimental resvilts that are available for the tribrid rocket motor indicate that the 
temperature variation with respect to z is small, and therefore, 02 t/ Bz^ can be 
eliminated. 

It was decided to use a backward difference numerical solution that is uncondition- 
ally stable (refs. 12 and 13) for the heat -transfer equation expressed in Cartesian coordi- 
nates coupled with a method that would insure convergence. From the computer runs it 
appears that no stability problems exist when the more complex cylindrical coordinates 
are used. The backward difference is similar in formulation to the explicit method 
except that the derivatives are calculated at the future time level. In order to do this 
calculation, equation (19) was put into finite -difference form and the temperature of a 
block at the present time t was expressed in terms of its temperature and the temper- 
ature of the surrotmding blocks at the future time t + At. 



^ ^ '^i.t + At - Tj,t 

at At 


(20a) 


a^T _ 1 fr^ 

ap2 Av 2\ (i+l)t + At 


8T _ 1 U 

ar 2 Ar\ (i+l)t + At 


“ ^'^i,t + At 




(20b) 

(20c) 
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Equations (20a), (20b), and (20c) can be substituted into equation (19) to obtain 


T,. 




r^2 Ar)(^(i+l)t+At " ’^(i-l)t+At^ 


(20d) 


and solving equation (20d) for yields 


* (2*i*Ar ■ t At + ^^2 


(*"^\t+At-(ff*^)T(W)t.At (2« 


Now a set of N simultaneous equations with N unknowns for each time step At can 
be formed. This equation is the general equation for the nonsteady conduction into a 
cylindrical solid. At this point the necessary boundary conditions that represent our 
particular problem must be imposed. 

Because the basic equation is parabolic, initial conditions, t = 0, and boundary 
conditions at r = Rq and r = R^g (fig. 3) must be known. The present problem also 
consists of defining the boundary conditions at the interface between the solid and liquid 
or at r = Rj^. The initial conditions (t = 0) will be assumed as: 

A very thin melt layer of at least four blocks. The first block will be at the 
surface temperature Tq; the next three blocks will decrease linearly from 
To to the melt temperature Tj^j. The blocks from the fourth block to the 
back-insulated wall R^g will be at some approximate ambient tempera- 
ture T^. 

The surface boundary condition is established as follows: 


Q 


(conducted into siirface) 



>1 


r=R 


o 


1a '^i.t + At " "^-l.t + At I 

2Ar ) 


( 22 ) 




J 


Because mass removal at the surface exists and must be taken into accoimt, an 
energy balance at the surface may be written 
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Qin = -k|^ + Pff hv 


^(radiation and convection to the surface) ^(conducted from surface) ^(vaporized) 


Qin “ Qcs ^ 


The term can be obtained from equation (17), and 


Qv = Pfrhy 

where 

Pj density of fuel 

r rate of surface regression with respect to time 

hy heat of gasification of fuel 

Then 

Qcs ~ ^in “ Pf^^v (23) 

Equation (23) can be substituted into equation (22) to obtain 

'^-l,t + At = '^l,t + At + (^in ■ Pf^*^v) (24) 

Equation (24) can be substituted into equation (21) and the resulting equation is the first 
equation in the set of simultaneous equations. The procedure is as follows: 

At the surface equation (21) takes the form 


AtK At 


“ \2ro Ar " Ar2 


+ 


. , 2 AtK^ 


+ 


AtK , AtK \ m /9-x 

l^j.2 2roArFl,t + At ^^5) 


Substitute equation (24) into equation (25) and let 
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AtK 


2rj) Ar 


b = 


AtK 

Ar2 



(a - b) t Qcs) 


+ (1 + 2b)To t - (a + b)Tj t + 


This equation reduces to 


To,t = (1 + 2b)To_t + 4t - 2»Ti_t + 4t + ' b)^Qcs 


(26) 


Because of the numerical method used, Ar at the surface must be much smaller than 
Fq and therefore b, the term with Ar2 in the denominator will be much greater 
than a. It therefore appears to be a safe assumption that the term a in equation (26) 
can be neglected. Equation (26) takes the form 


^o,t 


= 1 + 


2 At 
Ar2 




At 


2 AtK 


2 At 


^r2 l,t + At 


^ (Qin ■ Pf 


(27) 


From equation (27) it is seen that a new unknown r has been added to the set of simul- 
taneous equations and now there are N equations and N + 1 unknowns. Therefore, the 
assumption is made that the surface temperature Tq is known, and that it is constant 
for every time step or Tq ^ = Tq ^ Experiment and theory indicate that this 

assumption is reasonable. Equation (27) becomes 


where 


■^o,t(K2 - l) = Kj^Qin - P^hy) + *^3'^l,t+At 


Kl 


2 AtK 
Ark 


K2 

K3 



2 Atk\ 

Ar2 J 


2 AtK 


Ar'^ 


(28) 
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Equation (28) reduces to 



Kl 

Kg - 1 



.K2 - V 


( 29 ) 


where the left-hand side of equation (29) is considered to be known; the unknowns are 
t+At Equation (29) is the first equation in the set of N equations. 

The second equation in the set of N simultaneous equations is based on equa- 
tion (21): 


^l,t “ -^l^Oji + At + + 


where 


Ai =AtK — ^ + 


Bi = l + 


Ar^ 2r j Ar 

2 AtK 
Ar2 


Cl = -At K(— ^ 


2ri Ar ^r2. 


and Tp j considered to be known. Therefore, 


'^l,t " '^I'^o.t + At ” ®i'^l,t + At + ^l’^2,t + At (20) 

The third equation is 

T2,t = ^2Tl,t+At + ^Tg^t + At + ^2T3,t+At (31) 


where 

= 
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From here on the equations take on a repetitive form imtil the melt-solid interface is 
reached. This general equation, for the (i-1) block, is 


"^(i-l)! “ -^(i-l)'^(i-2)t + At + + (32) 


where 


A(i-l) + 2r(i_i) Ar 

C(i-i) = -AtKi 


\2r(i_i) Ar ^j,2j 

At this point the solid-liquid interface must be considered. Note the following figure: 



Solid-melt interface 


The fact that m is a moving interface with heat absorption can be taken into accoimt by 
performing the following heat balance: 


• • 

^(conducted into m) “ ^(heat of fusion) ^(conducted out of m) 


dTj 


- k. 


r=r 


m 


3Tp 
s ~dr 


„ 9r V 

r_r 


where the subscript I represents a liqxiid and the subscript s represents a solid. 
Thus, 


Qc,in - Qm + Qcs 


(33) 
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where 


Qc,in - -^(T(m-1) “ 






Equation (33) becomes 


kAm 

Ar 


(®’(m-l) “ "^m) “ Pf^m^m ■*" " "^(m+l)) 


(34) 


Because Ar is small, it can be assumed that Ajj^ ~ At this point it is noted 

that the energy balances performed on the surface and at the solid-liquid interface are 
not in cylindrical coordinates as is the general equation (21). Therefore, the area change 
between the surface and the interface must be taken into account. This change is deter- 
mined by setting Agm-f^ce = 1 3,nd since 6 and z are constant for all r 


A 


m 



By taking this fact into account, equation (34) is solved for Tm- Here the solution is for 
Tm instead of T(ni+1); this procedure simplifies the equation 


"^m,! + At “ 2 ^(m-l)t + At ■*" 


fAvVoPh^^ 1 ^ 
y 2krm + 2 (»i+l)t + At 


(35) 


Now substitute equation (35) into all equations containing Tj^^ as follows, and thereby 
introduce the desired quantity fjjj, the regression rate of the solid-liquid interface. At 
this point the assumption is made that Tm is known and that it is a constant. 

By following this procedure for the (m-1) block. 


"^(m-l)! " A(m-l)T(ni-2)t + At + ^ + ^2" ^^'^(m-l)t + 

^ /^rroPhii^j^ C(m-l) 

/ ™ 2 ^(m+l)t+/ 


At 


(36) 
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For the mth block, 


'^m,t - + At + ■*■ {f T(m+l)t + At 

For the (m+1) block 

m. . _-^(m+l) T,. .A 

T(m+l)t 2 ^ ^(m-l)t + At + ^(m+1)^ J ^ 

+ + At + ^(m+l)t + At (38) 

Equations (29), (30), (38), and (37) must be included in the melt layer, and therefore, four 
blocks are necessary. The rest of the equations in the interior of the solid are of the 
general form of equation (21). 

Now consider the insulated back surface (designated by bs). Note the following 
sketch: 



The boimdary condition for an insulated back is that no heat is transferred, or 

Q--kfE| 

^ 'r-Rbs 


P^bs+1 " '^bs 


Ar 


= 0 


Therefore 


"^bs+l - "^bs 

To eliminate the temperature Ti^g+i from the system of equations, let 
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"^bSjt " ■^bsTbs-1 + (b + C^^s^Tijs 


(39) 


At this point the system of simultaneous equations are complete and they take the fol- 
lowing form: 




Kl A 

" ’’o.V K2 - 1 


®i'^l,t + At + *'l'^2,t + At 

= Tl,t-AlTo,t 


A2Tl.t+At + %T2,t + At + C2T3,t + At 

= T2,t 


AlT(i-l)t + At + +At + CiT(i+l)t + At 

= Ti,t 


-3)t+At + BiT(m-2)t + At + C(m-2)T(m-l)t + At 

1 

II 


A(m-l)T(m-2)t + At+ (Bi + -i^)T(m-l)t + At -C(„.l)^ 2kr„ T(„^.i)t +At 



(ato + -^)T(m-l)t + At - Bl(^ 2kr„ )‘‘™ (2 '^m)T(nH-l)t + At 

II 

3 


” "^(m+l)t 


A(m+2)T(m+l)t + At + BiT(m+2)t + At + C(m+2)T(m+3)t +At 

“ '^(m+2)t 


ANT(N-l)t+At + BiTN,t+At + CNT(N+l)t+At 



where 


Abs-lT{bs-2)t + At + BiT(bs-l)t + At + C(bB-l)Tbs,t + At = T(bs-l),t 

(Abs - Cb^T(bs-l)t + At+ (Si + 2Cbs)Tbs,t + At = ^bs.t 


Kl 


2 AtK 
Ark 


K 2 = 1 + 


2 AtK 
Ar2 
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K3 = 


2 AtK 

Ar^ 


Ai = At K 


2ri Ar 



Bi = K2 


Ci = 



It is readily recognized that the system of equations could be put into a tridiagonal 
form were it not for the two equations preceding and following the equation for the melt 
layer. If the equations are modified and put into a tridiagonal form, their solution can be 
greatly simplified by using a recursion solution. Consider the equations (m-1), (m), 
and (m+1) simplified in the following manner for clarity: 


^l(m-l) + ^2(m-l) + C3(m-1) + C4(m-1) = 

® + ^2(m) + ^3(m) + ^4(m) +0 = 

^2(m+l) ^3(m+l) + C4(m+1) + C5(ni+1) = 


The terms that must be eliminated are and C2(m+1)' In order to eliminate 

the off-diagonal term C4(rn_i), each term in equation (m) can be multipled by 
C4(m-l)y/C4(m) and subtracted from each corresponding term in equation (m-1), with- 
out changing the value of the matrix. Equation (m-1) will become 


^l(m-l) " 


[ 0 ] 


,C4(m-l) 

^4(m) 


+ 


C2(m-1) - C2(m) 


^4(m-l) 
1^4 (m) 


^3(m-l) " ^3(m) 


^4(m-l) 

^4(m)_ 


+ 


^4(m-l) 

C4(m-1) - C4(m) 


- I^(m-l) " ^(m) 


^4(m-l) 

C4(m) 


which reduces to 


C2(m-1) - C2(m) 


C4(m-1) 

C4(m) 


+ 


C3(m-1) - C3(m) 


C4(m-1) 

^4(m) 


= °(m-l) 


■°(m) 


I^4(m-1) 
1^4 (m) 
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The same procedure is applied to equation (m+ 1 ), and it becomes 


C 3 (m+ 1 ) - C3(m) 


C 4 (m+ 1 ) ” C 4 (m) 


C2(m+1) 

C2(m) 


= D(m+1) 


-D(m) 


C2(m+1) 

^2(m) 


Now the set of simultaneous equations can be put into a tridiagonal matrix and solved by 
a simple recursion method. Appendix A gives the details of this procedure. 

In the matrix form the equations will be 


Cii + Ci 2 + Ci 3 

r 


i>o,t 

C22 + C23 + C24 

Tl,t + At 


Tl,t 

C33 + ^34 + ^35 

^m 

_ 

"^mjt 

C44 + C45 + €40 

« 

’^(m+l)t + At 

• 


'^(m+l)t 

• 

• 

• 

Ci,i + ^i,j+l ^i,i+2 

• 

• 

'^bs,t + At 


• 

• 

Tbs,t 


The solution of this set of equations will yield f, f and a temperature profile. From 
r and the thickness of the melt layer dm determined as follows: 

^m,t + At “ ^m,t + ^^*m “ 

From this equation, a time history of dm can be plotted. 

Determination of Physical Parameters 
Pertaining to Surface Heat Input 

In order to utilize a computer program developed by Marxman and Wooldridge at 
the Stanford Research Institute, the heat input to the fuel surface, certain parameters, and 
physical properties must be known. The most important of these parameters are also the 
most difficult to obtain; and a short discussion on how they were calcvilated is in order. 

It happens that the mass transfer number B incorporates most of the needed param- 
eters; B was previously defined as 
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where, from reference 2, 


and 


ue 



hv 


“ ^p,b^h “ ^p,f"^f,w 





where 

Ug 

— ratio of velocity in main stream to that at flame 

% 

AH sensible enthalpy difference between flame and fuel surface 

— oxidizer to fuel ratio 
F 

Ko, g concentration of oxidizer in free stream 

K mass fraction of metal or nonvaporizing material 

Cp^jj specific heat of gaseous products at reaction zone 

Cp^f specific heat of gaseous fuel at wall 

Tij flame temperature 

Tf,w temperature of fuel surface 

Tf Q ambient temperature of grain 

These parameters were obtained by calculations performed at Langley Research 
Center, from the manxifacturers of lithium and fluorine, and from references 14 and 15. 
The surface temperature was fovmd by experiment to be approximately 1700° F (926° C). 
(See ref. 13.) 



Physical Properties of the Solid Fuel 
The pertinent physical parameters of the solid fuel are 
k thermal conductivity state values 

p density 

Cp specific heat 

hy heat of vaporization 

hjjj heat of fusion 

All the preceding properties, except the thermal conductivity, can easily be obtained from 
handbooks, tables, and so forth. A method for obtaining the thermal conductivity of a 
nonhomogeneous substance such as the solid fuel is presented in reference 16. Refer- 
ence 16 deals with the effects of the addition of metal particles on the thermal conductiv- 
ity of plastics. In the solid fuel in question, metal particles are added to a hydrocarbon 
binder; thus, the problems are almost identical. Reference 16 indicates that two-phase 
systems consisting of powder and gas may be used to give preliminary thermal conduc- 
tivity data for two-component systems containing plastic and metal particles. The results 
were given in graphs where the abscissa is the ratio of solid conductivity to gas conduc- 
tivity and the ordinate is the ratio of effective conductivity to gas conductivity. A third 
parameter is the fraction of the space occupied by the gas. The thermal conductivity of 
the solid fuel system can be obtained by making the lithium metal analogous to the solid 
and letting the hydrocarbon binder be represented by the gas. The fraction of the space 
occupied by the hydrocarbon can easily be determined. By using this information in the 
graphs of reference 16, the thermal conductivity of the solid tribrid rocket fuel can be 
estimated. 


Computer Program Analysis 

The repetitive solution of a set of 200 to 300 simultaneous equations can prove to be 
rather time consuming when done manually. For this reason, a computer program was 
written whereby this same solution can be obtained in a very few minutes. The program 
is exceptionally easy to use and its use is a matter of suppl 3 dng a small number of input 
statements. Typical printouts include: temperature profiles through the grain; the 
regression rate of the surface RDOTV; the regression rate of the solid-melt interface 
RDOTM; the thickness of the melt layer DM; and the thickness of the solid layer DS. 
Any other output can be printed out if the user so desires. 
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The original computer program was written so that additional elements could be 
added to the melt layer as it increased and elements subtracted from the solid grain as 
it decreased. When the program was executed by using this method, it was found that 
when an element was added or subtracted, the melt-solid regression rate (RDOTM) 
became imstable and produced an erratic and unsatisfactory time plot for RDOTM. This 
condition was due to the fact that temperature gradients were changed discontinuously at 
those times, and the energy balances used in satisfying the boimdary conditions were, in 
turn, made discontinuous. This problem was solved by simply keeping the number of 
elements in the melt layer and solid grain constant and satisfsdng the convergence cri- 
terion that At = Ar^ for every calculation. Under these conditions the program pro- 
duces the expected form of time plots of regression rates and temperature profiles. 

It appears that the convergence criterion that At should be equal to or less than 
Ar^ is more than adequate. Comparison of computer runs by using this criterion on 
At and runs using a greater At (which is an indication of whether the solution is 
approaching the exact solution) shows no change in the eighth decimal place. This accu- 
racy is not necessary for engineering purposes. This matter is discussed because the 
program is taking from 10 to 15 minutes to analyze a 3-inch-thick cylindrical grain. 
Although this amount of computer time can be tolerated, it could be decreased by an 
order of magnitude that would be significant. 

Appendix B gives the necessary input and a listing of the program and some typical 
output. 


PRESENTATION AND DISCUSSION OF ANALYTICAL RESULTS 

The resvilts of this investigation are presented for: 

(1) slab burner 

(2) cylindrical grain 

(3) hypothetical grain 

The convective and radiative heat input to the fuel surface with the various fuel composi- 
tions was calculated for the cylindrical grain by using a computer program developed by 
Stanford Research Institute. This information was used as an input for the computer 
program that was developed as a resvilt of this study. Calculations were made for small 
laboratory motors with burn times ranging from approximately 20 seconds for fuels with 
the low-percentage Li grains (70 percent) to 10 seconds for fuels with the high-percentage 
Li grains (90 percent). For the slab burner an average value of the cylindrical grain 
heating rates for the same fuel formulation was used. In the case of the hypothetical 
grain, the heat input that had been calculated for the cylindrical grain was reversed with 


26 



respect to time; that is, the heat input that the cylindrical grain was experiencing at the 
end of burn time was made the heat input to the hypothetical grain at the beginning of 
burn. The heating rates used in each circumstance, as determined by the Stanford 
Research Institute program, are shown on the pertinent figures. 

Slab Burner 

In a slab burner the grain is in rectangular form. It burns with little surface area 
change, and this situation approximates a condition of constant heat input to the grain 
surface. Figures 4 show the calculated regression rates of the surface f and the 
solid-melt interface r^ as functions of burning time. Note that when f is larger 
than r, the solid-melt interface is receding faster than the surface is vaporizing; there- 
fore, the melt layer will increase in thickness. Figure 4(a) represents a pure lithium 
grain. Note the relatively large difference between r and and the fact that in the 

burning time r and f do not become equal; equal rates would indicate an equilib- 
rium condition. This situation can be improved by mixing the lithium with a low thermal 
conductivity hydrocarbon fuel, and thus allow less heat to be transferred into the grain to 
produce the melt layer. When the amount of lithimn in the fuel is decreased, the amovint 
of heat generated in the reacting boundary layer is also decreased. Figures 4(b) and 4(c) 
show the effects that the addition of 10- and 20 -percent hydrocarbon binder, respectively, 
have on the regression-rate curves. It is apparent from these figures that f can be 
decreased and equilibrium conditions obtained faster by the addition of a binder. 

Figure 5 presents melt-layer thickness as a fimction of burn time for an 80-percent 
grain burned at constant heat input to the surface or conditions typical of a slab burner. 
Note that the melt-layer thickness approaches a constant value in an asymptotic manner. 
Previous studies by United Technology Center on the other fuel formulations indicate that 
stripping will occur at a melt thickness of aroimd 0.125 inch (3.17 mm) for laboratory 
size motors. Figure 6 shows maximum melt-layer thickness as a function of percent 
lithium for a 3-inch-thick (7.62 cm) grain. Figure 6 gives a percent lithium value of 
approximately 79 percent at a melt-layer thickness of 0.125 inch (3.17 mm). 

Cylindrical Grain 

The analysis was extended to a small cylindrical grain rocket motor with an internal 
diameter of 3 inches (7.62 cm) and an outside diameter of 4.5 inches (11.43 cm). Four 
fuel formulations were analyzed containing 90, 85, 80, and 70 percent by weight of lithium. 
The results are shown in figures 7 to 10. 

Figures 7 show regression rates for the different fuel formulations as functions of 
burn time. These curves differ in two very significant ways from the preceding regres- 
sion rate curves for the slab burner. First, both r and r drop rapidly with time. 
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Second, the regression rates do not tend to approach each other and therefore never 
reach an equilibrium condition. Both of these phenomena can be explained by the fact 
that the heat input to the fuel surface is a function of the inverse of the port area of the 
motor which increases because of a receding fuel surface and a condition of constant 
oxidizer flow. The surface regression-rate drop is a direct consequence of decreased 
heat input per imit area. Because it takes much less energy to melt the fuel than to 
vaporize it, this condition of decreasing heat input per unit area will have less effect, 
proportionally, on the melt-solid regression rate than on the surface regression rate; 
therefore, they will tend not to converge within the reasonable application of cylindrical 
grains. 

Figure 8, which was generated by using surface heat input that was decreasing with 
time, shows that the melt-layer thickness continually increases with time, and never 
reaches an equilibrium condition. This condition is the condition experienced in a cylin- 
drical grain at a constant oxidizer flow rate. The melt layer continues to build up at a 
fairly constant rate and this buildup makes the cylindrical grain most impractical, 
especially in larger motors that have greater burning time. 

Figure 9 is a plot of melt layer at the end of burning as a function of percent lithium 
in the fuel. As can be seen, even 70 -percent lithium has a melt layer that is well above 
the suggested value of 0.125 inch (3.17 mm). 

Figure 10 gives two sets of siuface and melt-solid interface regression rates. The 
reason that two sets of curves were generated is as follows: After the melt layer builds 
up to a critical depth on the fuel surface, molten fuel droplets will be pulled into the gas 
stream by viscous forces. The fuel that leaves the surface in this manner will not have 
absorbed the energy that it would have, had it vaporized. Curves B show the regression 
rates obtained by assuming that the fuel is vaporized at the surface, and that the heat of 
vaporization is that energy needed to change the liquid fuel to a gas. Curves A repre- 
sent the condition of all the molten lithium being injected into the gas stream. Under 
these conditions an effective heat of gasification of the fuel was calculated where the 
hydrocarbon binder was decomposed and the lithium was only brought to the surface tem- 
perature as a liquid. These two curves represent the limiting condition of the process, 
and the regression rates should be bracketed by them. 

The preceding resvdts indicate that the cylindrical grain will produce more melt 
layer than the slab burner. High-speed photographs of the slab-burner firings show the 
occurrence of stripping; therefore, it is apparent that the cylindrical grain will produce 
even greater melt-layer removal. However, the effectiveness of a downstream turbu- 
lator will determine to a large extent the efficiency obtained by a cylindrical motor. 
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Also, it should be noted that because of an insulated outside boundary condition, the 
calculations show a melt-layer buildup that proceeds more rapidly toward the end of 
burning. This condition is due to heat being built up at the insulated wall, and energy 
that had previously been used to increase the grain temperature will at this point be used 
to melt the grain. The situation also presented itself in some of the slab burnings, and 
was noticed in the high-speed photographs that were taken. 

Hypothetical Grain 

After studying these results, it appears that if the heat input to the surface can be 
increased as a f\mction of time through manipulation of oxidizer flow rate and grain 
geometry, the preceding adverse trends could be reversed. In order to verify this theory 
analytically, computer rims were made where the heat input to the surface was increased 
with burn time; this procedure is the reverse of the former rims made with the cylindri- 
cal grains. It should be made clear that this situation is purely hypothetical. Figures 11 
to 13 show the effect that increasing heat input with burn time has on the regression rates. 
Melt-layer thickness as a function of percent lithium in fuel is shown in figure 12. This 
plot shows a definite improvement over figure 13. Figure 13 is a plot of melt-layer 
thickness as a function of burn time for the hypothetical case where the heat input to the 
surface is increasing with time. This graph shows that the melt layer decreases as the 
heat transfer increases beyond a certain point. At the point where the r curve 
becomes greater than the r^^j curve, the melt -layer thickness begins to decrease. It 
may well be possible to optimize fuel grain design and gas flow characteristics in order 
to minimize melt-layer buildup. This effect could be of even more importance in large 
motors where the greater burning time would allow an equilibrium condition to develop. 

Influence of Surface Heat Input on Melt-Layer Development 

Some important deductions concerning the burning characteristics of a high per- 
centage metal hybrid or tribrid solid fuel can be made from the preceding results. So 
that a direct comparison of the effect of heat input on the melt layer can be made, a 
graph of melt -layer thickness as a function of time for the 80 -percent lithium fuel for all 
heat input conditions considered is plotted in figure 14. After analyzing figure 14, the 
assumption may be made that the melt-layer thickness may be minimized by maximizing 
the surface heat input, with the reservation that this assumption applies to equilibrium 
conditions. 
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CONCLUDING REMARKS 


This investigation was ;andertaken to delineate the effects that an important param- 
eter, fuel thermal conductivity, has on melt-layer buildup and regression rate. An 
analysis of the results indicates, however, many of the effects of the heat input to the 
fuel surface from boundary-layer convection and radiation. The original objective was 
accomplished by producing an analytical method that is dependent upon experimental 
results for supplying the fuel surface temperature. 

The numerical method and mathematical procedure used in this study appear to 
have performed very satisfactorily. The solution was stable in all computer runs made. 
The convergence criteria proved to be more than adequate when solutions obtained by 
using different time increments were compared. 

A parabolic temperature profile is obtained rapidly, after a linear temperature 
profile is initially given to the melt layer and solid grain. The experimental results on 
regression rates available from slab burner firings agree with the calculations excel- 
lently. Motion pictures of the slab-burner firings show melt-layer buildup and the 
stripping phenomenon as predicted by the calculations. 

As expected, the results show a substantial decrease in the melt-layer thickness as 
the amount of hydrocarbon binder is increased in the tribrid fuel. Also, it was indicated 
that the surface heat input substantially affects the melt -layer buildup. It is recommended 
that a study be done to determine to what extent the surface heat input can be controlled in 
order to minimize the adverse effects of a thick melt layer. The surface heat input is a 
function of grain geometry and mass flow rates, these among other parameters could pos- 
sibly be constructively manipulated. 

Even though the physical properties for the fuel grain were obtained by using rea- 
sonable techniques and assumptions, experimental results would improve the overall 
accuracy of the method. This statement is especially true of the important parameter 
thermal conductivity. 

It is recommended that this analysis be applied to the problem of determining the 
aerodynamic and radiative heat input to the fuel surface. Using experimentally deter- 
mined regression rates and solving one set of equations for heat input would permit this 
application and would appear to be a possible approach for establishing the validity of our 
present theories for calculating surface heat input in many types of solid fuel motors. 

Langley Research Center, 

National Aeronautics and Space Administration, 

Langley Station, Hampton, Va., December 2, 1969. 
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APPENDIX A 


SOLUTION OF TRIDIAGONAL SET OF SIMULTANEOUS EQUATIONS 
The equations of interest are each of the form 


AjjTjj_j + = D^ (k = 1>2,. . . K) (Al) 

and constitute K equations in the K miknowns Ti,T2,. This set of equations 

may be represented by the matrix equation 

MT = D (A2) 

where T is the vector of the unknowns, M is the matrix operator, 



Pbi 

Cl 

0 

0 


A2 

B2 

C2 

0 


0 

As 

BS 

Cs 

M s 

0 

0 

A4 

B4 


and D denotes the resulting matrix. 

The matrix M is now factored i 


0 

• • • 

• 

• 

• 

0 

• • • 

• 

• 

• 

0 

« • • 

• 

• 

• 

C4 

• • • 

• 

• 

• 

• 

• « » 

• 

• 

• 

• 

• • • 

%-l 

Br-1 

Cr-I 

. 

• • 


■^K 

% 


two factors L and U; that is, 


(A3) 


M = LU 


(A4) 


which, because of the tridiagonal form of equation (A3) , may be taken in a particvilarly 
convenient form of a lower triangular matrix 
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wi 


U2 wg 0 


L s 


0 U3 Wg 


(A5) 


u 


K-1 ^K-1 


and an upper triangular matrix with unit diagonal 


Uk wk 


U = 


1 vj 0 


1 V2 0 


0 0 1 V3 0 


0 1 Vk_2 0 


'K-1 


(A6) 


Formal execution of the operations indicated by equation (A4) leads to the recursion 
formulas 


wi = Bi 


(A7) 


vi = 


wi Bj 


(A8) 


= Bk - AkVk_i 


(k = 2,3,. ..K) (A9) 


Ck 

Vk = — 
K Wk 


(k = 2,3,. ..K) (AlO) 
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UjjHAj^ (AH) 

so that the elements of L and U are readily determined. 

By definition, 


g= UT 


(A12) 


thus, eqviations (A2) and (A4) give 


LUT = Lg = D 


(A13) 


Then, by utilizing the lower triangular form of L, 



■ ^k^k-l 
Sk- Wk 


(k = 2,3,. ..K) (A14) 


By referring to equation (A12), only the elements of T yet remain and are given 
by the recursion formula 


Tk = Sk 


(A15) 


Tk = gk - VfcTk+1 (k = K,K-1,. . . 3,2) (A16) 


Thus, the unknown temperatures are generated from these equations in the order 

Tk»Tk- 1>- • -Ti- 

When the elements of the tridiagonal coefficient matrix depend upon the unknown 
temperatures (as is the case when properties are temperature dependent), an iteration 
procedure must be adopted in which the most recent value of the calculated temperatures 
is used in the computation of the elements of the matrix. The solution outlined is 
repeated successively until the computed temperatures agree within a desired precision. 
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APPENDIX B 


COMPUTER PROGRAM 

The necessary input and some appropriate comments are as follows: 

RO inside diameter of a hollow cylindrical grain (inches) 

BIGR outside diameter of a hollow cylindrical grain (inches) 

DMZ initial thickness of the melt layer. This value should be very small, on the 

order of 0.01 or 0.001 (inches) (0.254 mm or 0.0254 mm) 

DELRM thickness of a finite element in the melt layer. Because of the nature of the 
equations, there must be at least four elements in the melt layer (inches) 

DELRS thickness of a finite element in the solid grain (inches) 

DELT time increment between calculations. Reference 12 recommends a conver- 
gence criterion of At = AX^, This criterion appears to be more than 
adeqviate (seconds) 

TSURF surface temperature of grain (°F) 

TMELT melt temperature of grain (°F) 

HV heat of vaporization of the fuel (Btu/lb) 

HM heat of fusion Q of the fuel (Btu/lb) 

RHO density of fuel (ib/in^) 

CP specific heat of fuel (Btu/lb-^F) 

RT thermal conductivity of fuel (Btu/in-sec-°F) 

Q heat input to fuel surface due to convection and radiation (Btu/in^-sec) 

MO initial number of elements in melt layer 

TPRINT time increment on Data Printout (seconds) 
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TCR time at which calculations are to stop (second) 

TEMP(l) an initial temperature profile for the grain. The initial profile form is arbi- 
trary in that the proper configuration is rapidly assumed. This array is 
one -dimensional starting with the first value as the surface temperature 
and progressing to the back surface temperature 

The program listing is as follows: 


program D28S0 ( I NPUT « OUTPUT « T APE5= INPUT » T APE6 = OUTPUT ) 

REAL KD»KT»K1 «K2»K3 
INTEGER S 

D I MENS I ON R(300 )»A(300 ) »C(300> » DI AG 1 ( 300 > » D I AG2 ( 300 ) * DI AG3 < 300 ) « 
1 EL ! M ( 300 ) , D ( 300 ) » TEMP ( 300 ) * I DENT ( 8 ) 

2«G(300 ) «W(300) *V(300 ) tSAVET (300) 

DIMENSION QTAB ( 100 ) fROTAB < 100 ) 

EQU I valence (R(1 )«W(1 )).(A(1),V(1>)»<C<1)«G(1 ) ) 

NAMEL I ST/ I NPUT/RHO» HV » HM »KT * CP «RO « TMELT « TEMP « TSURF ♦ B I GR « DMZ 
1 « DELRM < DELT « DELRS »MO t KO » TPR I NT t Q i TCR «NQR » QT AB » ROTAB 
3000 READ(5*700 ) IDENT 

IF (EOF *5 ) 201 . 300 
201 STOP 

700 FORMAT (8A10) 

300 WRITE (6*701 > IDENT 

701 FORMAT ( IHl *8A1 0// ) 

READ (5* INPUT ) 

WRITE (6* INPUT) 

C 

C COMPUTE initial CONDITIONS 

ICHECK=0 
IERR=0 

KD=KT/(RHO*CP) 

DM=DMZ 
T I ME = 0 .0 
CHANM=0.0 
CHANS=0.0 
M=MO 
XM = M 

XT I ME = T I ME 

DS= (BIGR-RO)-DM 

K= (DM/DELRM )+ (DS/DELRS ) 

S = K-M 
RM=RO+DM 
GO TO 1000 
C 

C MAIN PROGRAM LOOP 

2000 RO=RO+RDOTV*DELT 

CALL FTLUP < RO * Q * 1 * NOR * ROT AB *QT AB ) 

TIME=TIME+DELT 
IF(TIME«GT«TCR)GO TO 3000 
DM=DM+ (RDOTM-RDOTV)*DELT 
rM=RM+ ( RDOTM*DElT ) 

IF (RO.GT.BIGR)GO TO 200 
C 

C COMPUTE NEW DELRM 
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DO 45 1=1 ♦< 

45 SAVET( I )=TEMP( n 
QUANT= (RDOTM-RDOTV)*DELT 

chanm=chanm+ouant 
IF (CHANM.LT.IOOO. ) GO TO 50 

C CHANGE THE NO. OF BOLCKS IN MELT LAYER WHEN THE THICKNESS HAS CHANGED .01 
C ADD 1 BLOCK for EVERY .01 ADDITION 
C GET NEW temperature PROFILE IN MELT LAYER 
MOLD=M 

CHANMA=CHANM/.05 
NEWBLK=CHANMA 
XM=YM+CHANMA 
M = XM 

CHANM=0 .0 

anew=newblk+i 

DIFF = (TEMP (MOLD-1 ) - TMELT) /ANEW 
DO 46 I=1.NEWBLK 
AI = I 

46 TEMP(M-I) = TMElT+A 1*DIFF 
K = M+S 

MPl =M+1 
I I =0 

DO 47 I =MP1 *K 
11=11+1 

47 TEMPI I )=SAVET(M0LD+I I ) 

ICHECK= ICHECK+1 

IF ( ICHECK.GT. 1 n )G0 TO 50 

144 WRITE (6,145) ( SA VET ( I ) , I = 1 , K ) 

WRITE (6,147) M, MOLD, NEWBLK,CHANl«IA,CHANM, DM, RM,RDOTM,RDOTV, TMELT 
1 ,K,DS,S , chans 

WRITE (6,146) (TEMP( I ) , 1=1 ,K) 

145 FORMAT (*0SA''ET*/(7tl8.7) ) 

146 FORMAT ( *OTEMP*/ ( 7E 1 8 . 4 ) ) 

147 FORMAT ( *0M = * I 5 , 3X , *MOLD = * I 5 , 3X , *NEWBLO = * I 5 , 3X , *CHANMA = *E1 4 . 6 , 3X 
1*CHANM=*E14.6,3X,*DM=#E14.6/ * RM=*E 1 4 . 6 , 3X , *RDOTM=»E 1 4 . 6, 3X*RD0TV 
2=+E14.6,3X,*TMElT=*E14.6,3X,*K= I*I 5,3X,*DS=*E 14.6,3X,*S=*I5/* CHAn 
3S=*ei4.6) 

50 DELRM=DM/FL0AT (M-1 ) 

C 

C COMPUTE DELRS 

C 

DS=BIGR-RM 

IF ((DS.LT. .1) • OR. (DM.LT. .005)) GO TO 8000 
X=RD0TM*DELT 
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CHANS=CHANS+X 

IF (CHANS. LT.IOOO.IGO TO 80 

C CHANGE NO. OF BLOCKS IN SOLID WHEN THICKNESS HAS CHANGED .01 
C SUBTRACT 1 BLOCK FOR EVERY .01 CHANGE 
C REMOVE BLOCK (MELT-1) 

CHANSA=CHANS/.05 

NCHANSA=CHANSA 

S=S-NCHANSA 

CHANS=0*0 

K = M+S 

MPl =M+1 

IF ( 1CHECK.GT.3 ) GO TO 80 
WRITE (6.145) (SAVET( I ) . 1 = 1 .K) 

WRITE (6.146) M .MOLD . NEWBLK . CHANMA . CHANM . DM . RM .RDOTM . RDOTV . TMELT 
1 .K.DS.S. CHANS 

WRITE (6.147) (tEMP( I ) . 1=1 .K) 

80 DELRS=DS/FL0AT(S ) 

K = M+S 

IF (K-300) 90.85.85 

85 WRITE (6.86) M . S ♦ DELT . CHANM . DS . CHANS 
RO=BIGR+l .0 

86 FORMAT (*0NUMBER OF BLOCKS EXCEEDS DIMENSIONS*/* M = * I 5 . 3X*S = * I 5 . 
1 3X*DELT=*E15.6.3X*CHANM=*E15.6.3X*DS=*E15.6.3X*CHANS=*E15.6 ) 

GO TO 40 

90 DELT= DELRM**2 

IF (DELRM.GT.DELRS ) DELT= DELRS**2 
1 OOO DO 2 I = 1 .M 

R ( I ) =RO+FLOAT { I )*DELRM 

A ( I ) =DELT*KD* ( ( 1 . / ( 2 . *R ( I ) *DELRM ))-(!. /DELRM**2 ) ) 

2 C ( I ) =-DELT*KD* ( ( 1 ./ (2.*R( I )*DELRM ))+(!. /DELRM**2 ) ) 

IF(DS.LT. .001 )GO TO 5000 

MP=M+1 

DO 3 I=MP.K 

R ( I ) =RM+FLOAT ( 1 -M )*DELRS 

A ( I )=DELT*KD* ( ( 1 ./(2.*R ( I )*DELRS ) )- ( 1 ./DELRS**2 ) ) 

3 C ( 1) =-DELT*KD* { ( 1 ./ (2.*R< 1 )*DELRS ) )+ ( 1 ./DELRS**2 ) ) 

C 

C COMPUTE diagonals OF MATRIX 

C 

5000 K1 = (2.*DELT*KD )/ ( DELRM*KT ) 

K2=l .+(2 .*DELT*kD)/DELRM**2 
K3= (2.*DELT*KD ) /DELRM**2 
DI AG2 ( 1 )=- (RH0*HV*K1 )/ (K2-1 . ) 

DIAG3I1 )=K3/(K2-1 %) 

DIAGl (2 )=0. 
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DI AG2 (2 ) =K2 
DI AG3 (2 )=C ( 1 ) 

MM=M-2 
DO 4 1=3, MM 

DIAGl ( I )=A ( I-l > 

D I AG2 ( I ) =K2 
DIAG3( n=C ( I-l ) 

4 CONTINUE 

IFIDS.LT. .001 )GO TO 6000 
DIAGl <M-1 )=A (M-2 ) 

DIAG2(M-1 ) =K2+C (M-2 )/2. 

DIAG3(M-1 )=-C(M-2)*( ( DELRM*RO*RHO*HM ) / ( 2 .*KT*RM ) ) 

EXTRAR=C (M-2 )/2, 

DIAGl C M ) = A ( M- 1 ) + ( K2/2 . ) 

DIAG2(M)=-K2*( (DE:LRM#R0*RH0*HM )/'<2.*KT*RM ) ) 

DIAG3(M) = (K2/2. )+C(M-l ) 

B=1 .+(2.*DELT*KD)/DELRS**2 
EXTRAL=(A(M)/2. ) 

DIAGl (M+1 ) =-A (M )* ( (DELRM*R0*RH0*HM ) / (2.*KT*RM ) ) 

DI AG2 (M+1 ) = ( A (M )/2. )+K2 
DI AG3(M+1 >=C (M ) 

M2=M+2 
KMl =K-1 
DO 5 I=M2,KM1 
DIAGl (I ) =A ( 1 -1 ) 

DIAG2I I )=B 
DIAG3( I )=C( I-l ) 

5 CONTINUE 

6000 IFCDS.LT. .001)K=M-1 
DIAGl (K )=A (K-1 ) 

DIAG2(K)=B+C(K-1 ) 

D1AG3(K)=0.0 

C 

C COMPUTE COLUMN VECTOR 

C 

D ( 1 )=TSURF- (K1 / (K2-1 . ) ) *0 
D ( 2 )=TEMP (2 )-A ( 1 )*TSURF 
DO 7 1=3, K 
7 D ( I ) = TEMP ( 1 ) 

IF(DS.LT. .001 )G0 TO 7000 
D(M)=TMELT 

D(M-1 IsDIM-l )-D (M )*EXTRAR/DI AG3 (M ) 

D(M+1 )=D(M+1 )-D (M )-»EXTRAL/DI AGl (M) 

PUT matrix into TRI-DIAGONAl form by ELIMINATING OFF DIAGONAL TERMS 


38 



o o o 


APPENDIX B 


D I AG3(M-1 )=DI AG3(M-1 ) - ( DI AG2 ( M ) *EXTRAR ) /D I AG3 (M ) 

DIAG2(M-1 )=DIAG2<M-1 )-(DIAG1 (M ) *EXTRAR ) /D I AG3 ( M ) 

DIAGl (M+1 )=DIAG1 (M+1 ) - ( D I AG2 (M ) *EXTRAL> /D I AG 1 (M) 

DIAG2(M+1 )=DIAG2(M+1 ) - ( D I AG3 t M ) *EXTRAL) /D I AGl (M) 

TR I -DIAGONAL MATRIX SOLUTION 

7000 W( 1 )»DIAG2( 1 ) 

V( 1 )=DI AG3 ( 1 )/W( 1 ) 

G( n=D( 1 )/W( 1 ) 

DO 20 JJ=2«K 

W< JJ) = DIAG2( JJ)-DIAG1 (JJ)*V(JJ-1 ) 

V ( JJ ) =D I AG3 ( J J ) XW ( J J ) 

20 G ( JJ ) = ( D ( J J ) -D I AG 1 ( JJ ) *G ( JJ- 1 ) ) /W ( J J ) 

KM=K-1 

TEMP(K)=G(K) 

DO 30 JJ=1 *KM 

30 TEMP(K-JJ)=GO<-JJ)-V(K-JJ)*TEMP(K-JJ+l ) 

RDOTV=TEMP ( 1 ) 

IF (RDOTV.lt. 0.0 )RDOTV=0.0 
RDOTM = TEMP (M ) 

IF (TIME.EQ.O. )GO TO 2020 
IF( <TIME-XTIME),LT. TPRINTIGO TO 2030 
XTIME=TIME 
2020 CONTINUE 

40 WRITE (6.705) TIME 
705 FORMAT (*0T I ME»:*E 1 5. 8/*0TEMPERA TURES* ) 

WRITE (6. 702 ) ( TEMP ( J J ) . J J= 1 . K ) 

702 F0RMAT(5(5X.E15.8) ) 

WRITE (6.704 ) RM.DM. DS ♦ DELRS . RDOTV . RDOTM . DELRM ,RO 
704 FORMAT (*0 RM =*E 1 5 . 8 . 5X . * DM=*E 1 5 . 8 . 5X . * DS=*E 1 5 »8 .5X .♦ DE 

1LRS=*E15.8/'* RDOTV=*E15.8.5X.*RDOTM=*E15.8.5X.*DELRM=*E15.8.5X.*RO 
2 = »E1 5.8 > 

WRITE(6.81 0 ) 

810 FORMAT (/* D ARRAY*/) 

WRITE (6. 702) (D( JJ) . JJ=1 .K) 

IF ( lERR.NE.O ) GO TO 3000 
2030 CONTINUE 

GO TO 2000 
200 CONTINUE 

WRITE(6.703 ) 

703 FORMAT (23H DM GREATER THAN (R-RO)) 

GO TO 3000 

8000 IERR=1 

GO TO 40 
END 
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TIMF= 2. 500H347E+00 


temperatures 

2.0t588970E-01 1 . 1 82 97415E+0 3 7 . 6fc369910E+02 2 . 0 6 1 7 19 15 E-0 1 2 . 307 1 2A72E+02 

1 .1 1436173E+02 1 .006 39992 E +02 1 .C0017356E + 02 1 . C0000422E + 02 1 .00000008E+02 

l.OOOOOOOOE+02 1. OOOOOOOOE+02 l.COOOOOOOE+02 

PM= 3 .62418303E+00 0M= 1 . 09950680E-02 DS= 8.7581&923E-01 DELRS= 9.73129914E-02 

R,f)0TV= 2.06588970E-01 RDOTM= 2 . 0 8 171 9 15 E-0 1 DELRM= 3. 65502266E-03 R0= 3 .6 131 880 lE+00 

D ARRAY 

-3.63252066E+03 1 . 18397993E + 03 7.665901^‘7E+02 3 . 50000000E + 02 2. 30933439E+02 

1.11436078E+02 1. 005399R3E+02 1 .CC017356E+02 1 . 00000422 E+02 1 .00000008E+02 

l.OOOOOOOOE+02 1 .OOOOOOOOE+02 1 .COOOOOOOE+02 

TIMF= 2.60012092F+00 

TEMPERATURES 

2.03833209E-01 1 . 18297056F +0 3 7.f-6366397E+02 2 . 054431 OlE-01 2. 31071 840E+02 

1. 12158601E+02 1. 00611547E+02 1 .00020966E+02 1 . 00000544E+02 l.OOOOCOllE+02 

1 .OOOOOOOOE+02 l.OOOOOOOOE+02 ' l.COOOOCOOE+02 

RM= 3.64486507E+00 DM= 1 . 1 1 547 16 1 E-0 2 DS= 8 .55 1 349 30E-0 1 0ELRS= 9 . 50149922E-02 

RDQTV= 2 .03833209E-01 RD0TM= 2 . 0 6443 1 OlE-0 1 r)ELRM= 3. 71 323 869E-03 RC= 3 .6 3371035 E+00 

D ARRAY 

-3. 63727678E+03 1 . 1 8397634E+03 7.£6586636£+02 3 .50000000E+02 2.3 1291803E+02 

1 .12158499E+02 1 .0061 1537E+02 1 .C0020966E+02 1 . 00000644E+ 02 1.0000001 1 E+02 

l.OOOOOOOOE+02 l.OOOOOOOOE+02 1 .COOOOOOOE+02 

TIME= 2.70012589E+00 

TEMPERATURES 

2.01114651E-01 1 .18296688E+03 7 .66362799E+02 2 . 0275 216 lE-01 2. 3 1444 267E+02 

1 .12909520E + 02 1. 006 90800E + 02 1 . C002 5 232 E+02 1 . 0000069 3 E + 02 1 .00000016E+02 

l-OOOOOOOOE+02 1 .OOOOOOOOE+02 1 . OOOOOOOOE+0 2 

RM= 3.66527555E+00 DM= 1. 1 31 70 879E-02 DS= 8 . 34724448E-C1 0ELRS= 9 .2747 1609 E-02 

ROOTV= 2.01114651E-01 RD0TM= 2 .02752151E-01 DELRM= 3.77236262E-03 P.0= 3 . 6 5395 846E+00 

D ARRAY 

-3.64222134E+03 1 .18397265E+03 7.66583039E+02 3. 50000000E + 02 2.31664225E+02 

1.12909411E + 02 1. 00690788F+02 1.C0025231E + 02 1 .000006S8E+02 1 .000000 16E+02 

l.OOOOOOOOE+02 l.OOOOOOOOE+02 1 .COCOOOOOE+02 
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Figure 3.- Cross section of tribrid fuel grain. 
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(a) 70-percent-lithiuni— 30-percent-binder; heat input to surface; Q = 12.7 Btu/in^-sec (20.75 x 10*5 w/m2) at burn time = 0; 

Q = 3.85 Btu/in2-sec (6.29 x 10& w/m^) at burnout. 

Figure 7.- Regression rates as function of burn time. Cyiindrical grain. 
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Bum time, sec 

(d) 90-percent-lithium— 10-percent-binder; Q = 45 Btu/in^-sec (?3.53 x 1Q6 w/m^) 
Q = 8.5 Btu/in^-sec (l3.89 x lO^ w/m^) at burnout. 


at time = 0; 


Figure 7.- Concluded. 
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Figure 10.- Regression rates as function of burn time. Cylindrical grain, 80-percent lithium-20-percent-binder. 

Dashed lines indicate extrapolation. 


Regression rates, In/sec 


o o o 



W1 


Regression rates, cm/sec 



Burn time, sec 


(a) 70-percent-lithium— 30-percent-binder; heat input to surface increasing with time; Q = 6 Btu/in2-sec ( 9.8 x 10® w/m2) 
at time = 0; Q = 12.7 Btu/in^-sec (20.75 x 106 w/m2) at burnout. 

Figure 11.- Regression rates as functions of burn time. Hypothetical grain. 






Bum tine, sec 

(d) 90-percent-lithium— 10-percent-binder; Q = 20 Btu/in^-sec (32.68 x 10® w/m^) 
Q = 45 Btu/in2-sec (73.53 x io6 w/m^) at burnout. 


at time = 0; 


Figure 11.- Concluded. 
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